/*
  non-Newtonian fluid with pressure- and shear-dependent viscosity
  
  Viscosity is of the form
   
     mu_0(A+|Du|^2)^((r-2)/2).
   
*/
#include <math.h>
#include "papcel/nnstokes3d.h"



C3DNonNewtonianStokesProblem::C3DNonNewtonianStokesProblem()
	: C3DFlowProblem()
{}


C3DNonNewtonianStokesProblem::~C3DNonNewtonianStokesProblem()
{
  outflows.clear();
  outflow_p.clear();
}

void C3DNonNewtonianStokesProblem::afterReadConfig(const Json::Value &input_obj)
{
  bool   setPressureMeanValue;
  int    pressureMVsegment = 0;
  double pressureMV = 0;
  
  if (input_obj.isMember("pressure_meanvalue"))
  {
    setPressureMeanValue = true;
    pressureMVsegment = input_obj["pressure_meanvalue"]["segment"].asInt();
    pressureMV        = input_obj["pressure_meanvalue"]["value"].asDouble();
    log(W_PROG, "Pressure meanvalue on %s = %e\n", mesh->getPhysicalName(pressureMVsegment).c_str(), pressureMV);
  }
  else
  {
    setPressureMeanValue = false;
  }
  
  registerQuadratures(&Q3d11p, &Q2d6p, &Q1d4p);
  addVector("velocity", FE_LIN);
  addScalar("pressure", FE_LIN, setPressureMeanValue, pressureMVsegment, pressureMV);
  

//	if (!input_obj.get("outflow", 0).isNull())
//	{
//		size_t begin = 0;
//		size_t comma = 0;
//		string bname;
//
//		while (comma < string::npos)
//		{
//		  comma = params["outflow"].find_first_of(",", begin);
//		  bname = input_obj["outflow"].substr(begin, comma-begin);
//		  outflows[mesh->getPhysNameId(bname)] = true;
//		  begin = comma+1;
//
//		  if (params.find("outflow_pressure[" + bname + "]") != params.end())
//			outflow_p[mesh->getPhysNameId(bname)] = atof(params["outflow_pressure[" + bname + "]"].c_str());
//		}
//	}

	// parameter convection indicates whether convective term is to be included into the equations
	if (input_obj.get("convection", "no").asString().compare("yes") == 0)
	{
		conv = 1;
		log(W_PROG, "Using convective term.\n");
	}
	else
	{
		conv = 0;
		log(W_PROG, "Not using convective term.\n");
	}

	// model of viscosity
	PAR_mu0   = input_obj.get("parameter_mu0", 1).asDouble();
	PAR_A     = input_obj.get("parameter_A", 0).asDouble();
	PAR_alpha = input_obj.get("parameter_alpha", 0).asDouble();
	PAR_r     = input_obj.get("parameter_r", 2).asDouble();
	PAR_q     = 2/(2-PAR_r);
	log(W_PROG, "\nViscosity parameters:\n");
	log(W_PROG, "mu_0  = %e\nA     = %e\nalpha = %e\nr     = %e\nq     = %e\n\n", PAR_mu0, PAR_A, PAR_alpha, PAR_r, PAR_q);

	// coefficients of outflow boundary conditions
	PAR_sigma1 = input_obj.get("parameter_sigma1", 0).asDouble();
	PAR_sigma2 = input_obj.get("parameter_sigma2", 0).asDouble();
	log(W_PROG, "Outflow coefficients:\nsigma1 = %e\nsigma2 = %e\n\n", PAR_sigma1, PAR_sigma2);

//	((C3DFlowProblem)(*this)).afterReadConfig(input_obj);
}

void C3DNonNewtonianStokesProblem::afterPrepare()
{
	char e[1023];

	// start with initial guess different from zero
	if (!loaded_solution)
	{
		for (int i=0; i<q.size(); q(i) = q(i) + 1e-0, i++);
	}
	else
	{
	//    for (int i=0; i<q.size(); q(i) = q(i) - q(i).getVal()*0.9, i++);
	}


	if (PAR_r != 2)
	{
	  sprintf(e, "%e*(%e+gvelocity00^2+gvelocity11^2+gvelocity22^2+(gvelocity01+gvelocity10)^2*0.5+(gvelocity02+gvelocity20)^2*0.5+(gvelocity12+gvelocity21)^2*0.5)^(%e)", PAR_mu0, PAR_A, (PAR_r-2)/2);
	  addPPscalar("viscosity", e);
	}
	else
	{
	  sprintf(e, "%e", PAR_mu0);
	  addPPscalar("viscosity", e);
	}

	addPPscalar("normDv", "sqrt(gvelocity00^2+gvelocity11^2+gvelocity22^2+(gvelocity01+gvelocity10)^2*0.5+(gvelocity02+gvelocity20)^2*0.5+(gvelocity12+gvelocity21)^2*0.5)");
}

void C3DNonNewtonianStokesProblem::afterEvaluateVolumeComponents(int elem, int qp)
{
  if (PAR_r != 2)
  {
	CVar Du2 = U[0]->grad(0)*U[0]->grad(0)
		+U[1]->grad(1)*U[1]->grad(1)
		+U[2]->grad(2)*U[2]->grad(2)
		+(U[0]->grad(1)+U[1]->grad(0))*(U[0]->grad(1)+U[1]->grad(0))*0.5
		+(U[0]->grad(2)+U[2]->grad(0))*(U[0]->grad(2)+U[2]->grad(0))*0.5
		+(U[1]->grad(2)+U[2]->grad(1))*(U[1]->grad(2)+U[2]->grad(1))*0.5;
	visco  = pow(Du2 + pow(pow(U[3]->val*PAR_alpha,2)+1,-PAR_q/2) + PAR_A, (PAR_r-2)/2)*PAR_mu0;
  }
  else
  {
	visco = PAR_mu0;
  }
  diverg = U[0]->grad(0) + U[1]->grad(1) + U[2]->grad(2);
  
  delta = pow(U[0]->ip3d->getDJ()*6, 2./3.);
}

CVar C3DNonNewtonianStokesProblem::volumeIntegral(int comp, int elem, int qp)
{
  switch (comp)
  {
    case 0:
    case 1:
    case 2:
      if (tau>0)
      {
        return (U[comp]->val*3-U[comp]->o_val.getVal()*4+U[comp]->oo_val.getVal())/(tau*2)*phi
                +(
                     (U[comp]->grad(0)+U[0]->grad(comp))*gphi(0)
		    +(U[comp]->grad(1)+U[1]->grad(comp))*gphi(1)
		    +(U[comp]->grad(2)+U[2]->grad(comp))*gphi(2)
		 )*0.5*visco.getVal()
            +(U[0]->val*U[comp]->grad(0) + U[1]->val*U[comp]->grad(1) + U[2]->val*U[comp]->grad(2))*phi*conv
            -U[3]->val*gphi(comp);
      }
      else
      {
        return visco*(
                    (
                     (U[comp]->grad(0)+U[0]->grad(comp))*gphi(0)
		    +(U[comp]->grad(1)+U[1]->grad(comp))*gphi(1)
		    +(U[comp]->grad(2)+U[2]->grad(comp))*gphi(2)
		    )*0.5
		   )
            +(U[0]->val*U[comp]->grad(0) + U[1]->val*U[comp]->grad(1) + U[2]->val*U[comp]->grad(2))*phi*conv
            -U[3]->val*gphi(comp) + diverg*gphi(comp)*delta
            -phi*(comp==2?-1:0);
      }
      break;
    case 3:
      return phi*diverg + (U[3]->grad|gphi)*delta;
      break;

  }
  return 0;
}



CVar C3DNonNewtonianStokesProblem::edgeIntegral(int comp, int elem)
{
  return 0;
}

void C3DNonNewtonianStokesProblem::saveOutflow()
{ // save velocity fields along the 'outflow' segments
  for (map<int, bool>:: iterator it = outflows.begin(); it != outflows.end(); it++)
  {
    ofstream f((params["savebase"] + "_" + mesh->getPhysicalName(it->first) + ".dat").c_str());
    for (int e=0; e<mesh->getNfaces(); e++)
    {
      if (mesh->getFacePhysicalNo(e) == it->first)
      {
        for (int i=0; i<3; i++)
	{
	  int ind = mesh->getFaceNode(e, i);
	  f << mesh->getCoord(ind, 0).getVal() << " "
	    << mesh->getCoord(ind, 1).getVal() << " "
	    << mesh->getCoord(ind, 2).getVal() << " "
	    << q(ind).getVal() << " "
	    << q(mesh->getNnodes()+ind).getVal() << " "
	    << q(2*mesh->getNnodes()+ind).getVal() << endl;
	}
      }
    }
    f.close();
  }
}

void C3DNonNewtonianStokesProblem::savePressure()
{ // save pressure along the 'outside' segment
  for (map<int,string>::iterator it=mesh->beginPhysicalNames(); it != mesh->endPhysicalNames(); it++)
  {
    ofstream f((params["savebase"] + "_pressure_" + it->second + ".dat").c_str());
    for (int e=0; e<mesh->getNfaces(); e++)
    {
      if (mesh->getFacePhysicalNo(e) == it->first)
      {
        for (int i=0; i<3; i++)
	{
	  int ind = mesh->getFaceNode(e, i);
	  f << mesh->getCoord(ind, 0).getVal() << " "
	    << mesh->getCoord(ind, 1).getVal() << " "
	    << mesh->getCoord(ind, 2).getVal() << " "
	    << q(3*mesh->getNnodes()+((QuadMesh *)mesh)->getNoOfLnode(ind)).getVal() << endl;
	}
      }
    }
    f.close();
  }
}







